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SUMMARY 

This report presents two numerical methods to describe the unsteady flow 
field In the blade-to-blade plane of an axial fan rotor. These methods solve 
the compressible, time-dependent, Euler and the compressible, turbulent, time- 
dependent, Navier-Stokes conservation equations for mass, momentum, and energy. 
The Navier-Stokes equations are written in Favre-averaged form and are closed 
with an approximate two-equation turbulence model with low Reynolds number 
and compressibility effects Included. The unsteady aerodynamic component is 
obtained by superposing Inflow or outflow unsteadiness to the steady conditions 
through time-dependent boundary conditions. The integration in space Is 
performed by using a finite-volume scheme, and the Integration in time is per- 
formed by using k-stage Runge-Kutta schemes, k = 2,5. The numerical integra- 
tion algorithm allows the reduction of the computational cost of an unsteady 
simulation involving high-frequency disturbances in both CPU time and memory 
requirements. Less than 200 sec of CPU time are required to advance the Euler 
equations In a computational grid made up of about 2000 grid points during 
10 000 time steps on a CRAY Y-MP computer, with a required memory of less than 
0.3 megawords . 


INTRODUCTION 

The reduction of automotive fan noise is now the subject of vigorous 
research. A better understanding of fan noise generation calls for a proper 
description of the physical processes which connect unsteady flow interactions 
with fan blades to far field noise. The unsteady aerodynamic component of the 
fan flow field is linked to the resultant far field acoustic pattern via the 
blade unsteady aerodynamic response, the coupling with the duct, the propaga- 
tion in the duct, and the acoustic coupling to the far field noise. The 
present work is devoted to numerically simulating the unsteady flow field 
through the blades of an axial flow fan rotor by solving the Euler and Navier- 
Stokes equations . 


*Work funded by Space Act Agreement C-99066-G. 


In a classical description of the fan flow field, the actual three- 
dimensional flow is built up by superposing two two-dimensional flow solutions 
to yield a quasi-three-dimensional flow field. The first solution is a 
through-flow solution, which treats the circumferentially averaged flow. The 
through-flow solution defines the flow conditions upstream and downstream of 
the rotor blades on axi symmetric surfaces. The second solution is a blade- 
to-blade flow solution on these surfaces, defining the flow conditions Inside 
the rotor blade passage. 

In a design approach, simplified through-flow and blade-to-blade solutions 
plus empirical correlations are used to define the shape of the airfoils along 
the blade-to-blade surfaces. These airfoils are then staked to generate the 
blade. The staking includes blade forward skew resulting from simplified 
acoustic solutions for the Isolated airfoil unsteady aerodynamic response plus 
simplified solutions for the acoustic coupling to the far field. The resulting 
fan geometries are therefore tested experimentally in order to assess the over- 
all performances. 

In an improved design approach, analytical methods become useful to intro- 
duce changes in blade shape. These changes are usually defined by using steady 
flow analysis methods in order to improve steady blade performances. The use 
of unsteady flow analysis methods should be able to improve both steady and 
unsteady blade performances. The present report describes two computer codes 
developed to study the steady and unsteady flow on the blade-to-blade surfaces. 
The flow field is predicted by numerically solving the compressible, time- 
dependent, Euler and the compressible, turbulent, time-dependent, Navier- 
Stokes conservation equations for mass, momentum, and energy. The areas of 
major concern in solving the governing conservation equations are the modeling 
of turbulence in the Navier-Stokes conservation equations, the formulation and 
implementation of appropriate boundary conditions, and the formulation of a 
time-accurate numerical integration algorithm. 

The Navier-Stokes equations are written in Favre-averaged form. The con- 
servation equations for mass, momentum, and energy are closed by using a two- 
equation turbulence model with low Reynolds number and compressibility effects 
included. The turbulence model is valid only for steady flows, but despite the 
lack of theory, the prediction capability of two-equation turbulence models in 
unsteady flow computations is not worse than the one obtained in steady flow 
computations. The proposed turbulence model suggests wall and compressibility 
corrections to standard K-e turbulence models in order to extend their range 
of applicability to compressible and separated flow fields. Apart from the 
values of the model constants actually specified, the model should be able to 
properly predict the flow turbulence when experimental data are available in 
enough detail to guide the model development and fix the values of the con- 
stants better. 

The unsteady aerodynamic response is obtained by superposing inflow or 
outflow unsteadiness to the steady conditions. For an unexcited jet flow, the 
inflow total pressure and density and the outflow static pressure are specified 
at constant values. To facilitate investigating the influence of acoustic or 
pressure excitation, the inflow total pressure and total density or the out- 
flow static pressure are modified by imposing time-dependent values. The time 
dependence is usually accomplished by imposing a sinusoidal oscillation to the 
steady values or a simpler step function. For now, the unsteady blade 
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performances can be defined in terms of blade pressure fluctuation amplitude 
normalized by the forcing pressure perturbation amplitude, but further work is 
needed in order to better define these unsteady blade performances. 

The integration in space is performed by using a finite-volume scheme. 
Along the boundaries, all of the flow parameters not specified by a boundary 
condition are computed by solving their conservation equations over half of 
the control volumes. The integration in time is performed by using multistage 
Runge-Kutta schemes. The numerical integration algorithm has been developed 
with the aim of reducing the computational cost of an unsteady simulation 
Involving high-frequency disturbances. The finite-volume formulation is 
adopted to allow the use of an arbitrary grid in the application to complex 
flow geometries. The multistage Runge-Kutta scheme is used both to improve 
accuracy in time and to extend the stability region of the explicit discretiza- 
tion. Two-, three-, four-, or five-stage Runge-Kutta schemes similar to 
Jameson's (Jameson et al . , 1981) are used. A k-stage scheme can be made stable 
for a Courant number up to k - 1 , but these schemes cannot be more than 
second-order accurate in time (Chima, Turkel, and Schaffer, 1987). The order 
of the scheme can be selected in accordance with the stability requirements 
for the particular problem. Dissipative terms have been added to suppress the 
tendency for odd and even points decoupling and to prevent the appearance of 
wiggles in regions containing severe pressure gradients. The numerical dissi- 
pation is introduced as the difference between second- and fourth-order differ- 
ences. The magnitude of the dissipative terms has been adapted to the local 
flow properties by means of a sensor based on the local pressure gradient. 

The overall accuracy is second order in both space and time. 


EULER EQUATIONS 

The computation of high Reynolds number, unseparated, unsteady cascade 
flow fields can be easily accomplished by solving the Euler conservation equa- 
tions written in unsteady, compressible form. The basic variables are the gas 
density, velocity, and enthalpy. Their conservation equations are written as 
fol lows : 

p.t + <f> V.k ■ 0 

(p u,) it * (p u k U,> k - -P >t 
(p h 0 )_ t + (P Uj h 0 ) j - P it 

where the state equation is 

p = ( y - 1 ) /y p • h 

The integration of the previous equations requires the introduction of the 
conditions to be fulfilled along the boundaries of the domain of interest. The 
physical flow domain is represented in figure 1. (As an aid to the reader, a 
symbols list has been included in the appendix.) The flow domain Is limited by 
inflow, outflow, solid, and periodic boundaries, where appropriate boundary 
conditions are needed. Along the inflow boundary, the total pressure poi . 
total density poi , and flow angle 3i are specified, whereas the Mach number 
is extrapolated from the Interior. These values can be both steady or time 
dependent. Along the outflow boundary, the static pressure P 2 is specified, 
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whereas all other variables are extrapolated from the interior. The value 
used for the static pressure can be steady or time dependent. The inflow and 
outflow boundary-condition values are assumed to be constant in the pitchwise 
direction, as theoretically fulfilled only at the upstream and downstream infi- 
nite. Across the solid boundaries, the flux of mass and energy is set equal 
to zero, whereas the flux of momentum is evaluated from the blade pressure and 
all of the variables are extrapolated from the Interior. 

Before starting the time integration, proper initial conditions have to 
be specified. These conditions can be either a previously computed steady- 
state condition or a simple no-flow condition. 


NAVI ER-STOKES EQUATIONS 

A better description of the unsteady flow field can be accomplished by 
solving the Navi er-Stokes conservation equations. The solution of the Navler- 
Stokes equations is required in the study of the cascade flow field at off- 
design conditions, when the incidence effects play an important role. The 
full Navi er-Stokes conservation equations are written in unsteady, compressible 
form. The basic variables are the gas density, velocity, and enthalpy. Their 
conservation equations are written in Favre-averaged form as follows: 


p,t * ( p Wk ■ 0 


<p u, ) t ♦ <P u k V >k - -p_, ♦ (i tk - P u" k u”,) i| 


(p v >t * cp u 3 v.j - p.t * u "j p.3 * t-«j - p f'l.j 


* [ "u 3 t ’ p u "i u "j 5 i ! .j * u ”t,j T u * p u "i u "k 5 i,k 
where (Hirsch and Deconlnck, 1985) 


h 


0 “ 


h + 


u u 


i 1 
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p u " u " 
1 1 


2 P 


and the state equation is 


P = (y - l)/yp*h 

and the constitutive relations are 

* '“(“l.J + “3.1 - I 5 k,k d ij) 


and 
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The previous equations are not closed, and some empirical expressions 
have to be Introduced for the turbulent quantities. A turbulence model which 
reflects the availability of experimental data Is needed. The closure Is 
achieved by using a first-order turbulence model. These models are probably 
the most reliable turbulence models because of the availability of experimen- 
tal results and computer resources. The K-c turbulence model adopted here 
Includes low Reynolds number terms (to account for near wall effects and to 
compute separated flows, In accordance with the suggestions given by Borettl 
and Mar te 111, 1989a) and Includes a pressure-density gradient term (to account 
for compressibility. In accordance with the suggestions of Jones and McGuIrk, 
1987). Note that the model Is strictly valid only for steady flows but that 
it can be used In unsteady flow computations without Introducing great errors. 

The Reynolds stresses and the turbulent heat flux vector are expressed as 
fol lows : 


p u", ■u-j - -P t (u u * u j-t - f <f tj ii,,,) ♦ f d tJ Kp 


P u'Y h" - 


->*t h .t 


Pr, 


where 


p c f r 
u u 

^t — t _h 


The Favre-averaged conservation equations for K and e are written as 
fol lows : 


(p K) t + (p u 1 K) i 




+ P + F - p c + 0 


,1 


(p e) £ + (p U 1 c) | 


Sc g (p t + p) 


*■ H C el f , * 6> - C c2 f 2 p * E 


,1 


where the turbulence kinetic energy production Is 


P = -p u". u". u. . 
p i 1 i.j 

the wall terms D and E are (Borettl and Mar te 111, 1988) 

,2 
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and the compressibility terms F and G are (Jones and McGuirk, 1987) 



The compressibility terms F and G play an Important role only for 
high-speed flows, and they have a negligible influence at low speeds. 

The values of the model constants are chosen as follows: 

Cjj = 0.09 C e i = 1.43 C e 2 - 1-92 Sc< » 1.0 Sc e = 1.3 
and the low Reynolds number functions are expressed as 

fy = expC-2 .5/(1 + R t /50)] 
f] = 1.0 
and 


f 2 = l.o - 0.33 exp[-R t 2 ] 
where the turbulent Reynolds number is 


Final ly, 


and 



U = -p C 

The boundary conditions for the Navler-Stokes equations differ from those 
adopted for the Euler equations because of the Influence of viscosity. Along 
the Inflow boundary, the total pressure pqi , total density poi, flow angle 
3], Inlet turbulence level Tu i , 

K - Tu V 2 

and the length scale of the turbulence motions Lei, 
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_ k 3/2 

e ' Le 


are specified, whereas the Mach number is extrapolated from the interior. 

These values can be both steady or time-dependent values. Inflow boundary- 
condition values constant In the pitchwlse direction are theoretically ful- 
filled only at the upstream infinite, but the influence on the cascade flow is 
not too strong if enough mesh points are used upstream of the blade cascade. 
Along the outflow boundary, the static pressure P 2 is specified, but all of 
the other variables are extrapolated from the interior. The value used for 
the static pressure can be a steady or a time-dependent value. An outflow 
boundary-condition value constant in the pitchwlse direction is theoretically 
fulfilled only at the downsteam infinite but has no Influence on the cascade 
flow when enough mesh points are used downstream of the blade cascade. Along 
the solid boundaries, the velocity, turbulence kinetic energy, and turbulence 
kinetic energy dissipation rate vanish. For adiabatic flows, the temperature 
gradient normal to the wall is set equal to zero. Otherwise, the temperature 
is set equal to the wall value. Across the solid boundaries, the flux of mass 
is set equal to zero, along with the flux of energy for adiabatic flows. The 
flux of momentum, the flux of energy for nonadiabatic flows, and the fluxes of 
the turbulent variables are eventually evaluated from the blade pressure and 
viscous terms. The density is finally extrapolated from the Interior. 


NUMERICAL METHOD 

The computational domain is shown in figure 2. The grid is obtained by 
the intersection of pseudostreamlines and pitchwlse lines. The equations are 
discretized in space by using a cell-centered, finite-volume formulation. 

Figure 3 shows the finite-volume molecules approximating Inviscid and viscous 
terms. The time derivatives and the source terms are evaluated at the cell 
center grid point. The spatial differences are evaluated as net fluxes across 
the faces of the cell, computed at the cell faces, and divided by the cell 
area. The discretization in space of the inviscid terms of the Navi er-Stokes 
equations < 1 . e . , the Euler equations) thus requires the cell center grid point 
and the 6 surrounding grid points. Since the first-order derivatives appear- 
ing in the viscous terms of the Navi er-Stokes equations (stress tensor, heat 
flux, and others) are expressed by using 4 grid points, the discretization in 
space of the Navler-Stokes equations involves 17 grid points. Global conserva- 
tion is ensured by evaluating the flux vectors on the faces of the boundary 
conforming mesh cell. Half-cells are used along the inflow, outflow, and 
solid boundaries. The first-order derivatives appearing in the viscous terms 
of the Navi er-Stokes equations (stress tensor, heat flux, and others) and the 
first-order derivatives of the inviscid terms are expressed along the bounda- 
ries by using one-sided differences. The scheme is second-order accurate 
in space on a regular grid. 

The equations are Integrated in time by using explicit, multistage, dis- 
sipative, Runge-Kutta algorithms. When the previous time-dependent equations 
are rewritten in the following vector form 

f , = T(f ) 

» t 
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the solution at the time t = (m + l)dt Is expressed as a function of quanti- 
ties evaluated at the times (m)dt according to the following scheme: 


f 


( 1 ) 


^(m) 


^(0) = ^m 

= f (0) - dt [T i <f <0> ) + T v (f (0) ) + AV 2 D 2 (f (0) ) - AV 4 D 4 (f (0) )] 
f (0) - C k dt CT 1 (f <k “ 1> ) + T v (f (0> ) + AV 2 D 2 (f (0> ) - AV 4 D 4 (f <0) )] 


^(m+1 ) = f (m) 

where D 2 Is a term approximating second-order differences, D 4 Is a term 
approximating fourth-order differences, AV 2 Is a viscosity coefficient for 
second-order differences, where 

AV 2 = C ? abs(D 2 (p))/(2p) 

(AV > 

with p an averaged pressure, and AV 4 Is a viscosity coefficient for 
fourth-order differences 


AV 4 - max(0.0, C (AV 4 } - AV 2 ) 

The second-order dissipative term as well as the averaged pressure Is 
evaluated in every mesh cell by using the 7-grid-polnts scheme of figure 3. 
Since the fourth-order terms are expressed as second-order differences of 
second-order differences, the discretization In space of these terms Involves 
the 19 grid points. The residual vector is split into two parts, an Invlscid 
and a viscous one, for convergence-improving purposes. The artificial vis- 
cosity coefficients are chosen as follows for Navler-Stokes solutions: 


C (AV 2 ) 


0.27 


C (AV 4 ) 


0.016 


whereas slightly larger coefficients are used for Euler solutions. 

The time step limit for an explicit scheme is established by the classi- 
cal Courant-Frl edri ch-Levy stability limit (CFL) for Invlscid flows. 


dt < dr 


CFL 
a + V 


where dr is a characteristic dimension of the space discretization, a is 
the speed of sound. 


2 - - 
a = Y* p/p 

and V Is the velocity. The order of the multistage Runge-Kutta scheme Is 
selected In accordance with the time step and stability requirements. By 
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Increasing the order of the multistage Runge-Kutta scheme, the maximum CFL 
number increases, and a proper selection of the scheme order follows from the 
prescribed time scale of the unsteady flow. A k-stage scheme used with cen- 
tral differencing is stable for a Courant number of approximately k - 1, 
depending on the choice of Cj , 1 = l,k. The influence of a CFL greater than 

unity on the maximum time step is shown in figure 4 for a CFL of about 2. For 

consistency, C^ = 1 , whereas for second-order accuracy in time C|<_i = 1/2. 

Note that the schemes of this form cannot be more accurate than second-order 
for any values of Cj , i = l,k (Chima, Turkel, and Schaffer, 1987). The values 
of Cj, 1 = l,k proposed by Chima et al . , 1987 for two-, three-, four-, and 
five-stage schemes are adopted here, allowing CFL numbers of about 1.00, 1.50, 
2.75, and 3.50, respectively. A switch parameter enables the choice of the 

k value. A CFL smaller than the one permitted by a linear inviscid stability 

analysis is adopted in the solution of the Navier-Stokes equations. This 
choice was made because of the full coupling of the K-e conservation equations 
with the basic variables conservation equations and the neglected viscous 
effects contribution to the definition of the maximum time step and also 
because of the stronger mesh size variations over the computational domain. 


RESULTS 

The numerical method previously described was applied to the computation 
of both steady and unsteady flow fields. The steady computations were per- 
formed to assess the basic code steady-flow prediction capability, and the 
method was applied to the computation of the flow within the low-turning fan 
rotor blade section shown in figure 5. The blade was designed by J. Sanz 
(Sanz, 1988, 1989) for an inlet Mach number of Mj = 0.20 and an inlet flow 
angle of 3j = 70.00° (with reference to the axial direction). The blade pro- 
vides low turning, giving an outlet flow angle of 32 = 65.00° and reducing 
the flow speed to an outlet Mach number of M 2 = 0.16. The pi tch-to-chord 
ratio is 0.75, the maximum thickness-to chord ratio is 0.05, and the blade 
chord i s about 0.10 m. 

The steady-flow computation was performed as a transient computation, 
starting from Initial no-flow conditions and suddenly decreasing the outflow 
static pressure to the prescribed steady value. The required consistency in 
time reduces the convergence properties with respect to those obtained in 
pseudo-unsteady methods (Boretti, 1989b), but the aim of this analysis is the 
simulation of an unsteady flow behavior, not better convergence to a steady 
state. 

At the beginning of the calculations, the Inflow and outflow boundary con- 
ditions are given as follows 

PQj = 100 000 N/m^ 

Pq ] = 1.21 kg/m^ 

Bi = 70.00° 

Tuj = 0.05 
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Le] = 0.005 m 
P2 = P01 

The outflow static pressure Is then suddenly reduced to the steady value 

P 2 = 98 000 N/m 2 

and the flow starts to exit from the computational domain. After a certain 
number of time steps, a new steady-flow condition is reached. 

The unsteady flow was simulated on the NASA-Ames Cray Y-MP. The computa- 
tions were performed by using time steps of the order of 10~ 4 sec, close to 
the time step limit with CFL equal to unity. The two-stage scheme was 
selected during the computations, thus reducing the total computational 
effort. The advance of the flow equations in a computational grid made up of 
about 2000 grid points performing 10 000 time steps requires about 300 sec CPU 
time when solving the Navier-Stok.es equations and about 200 sec CPU time when 
solving the Euler equations. The memory required is less than 
0.3 megawords . 

. First, an Euler solution was performed. The computation mesh adopted in 
the inviscid computations is presented in figure 6. The grid is a 25 by 75 
H-grld, made up of 25 computational points in the pitchwise direction and 75 
computational points in the axial direction, wi.th 41 computational points 
between leading and trailing edges. Figure 7 shows the computed surface Mach 
number distributions. The predicted flow field is very close to that 
expected, with a large velocity peak close to the blade leading edge followed 
by a smooth flow deceleration on the suction side and an almost constant veloc- 
ity distribution on the pressure side. 

The flow deceleration on the blade suction side can be strong enough to 
produce an undesired flow separation. In order to check the flow behavior on 
the rear suction side, a more accurate Navier-Stokes solution has been 
performed. The computational mesh adopted in the viscous computations is pre- 
sented in figure 8. The grid is a 35 by 85 H-grid, made up of 35 computa- 
tional points in the pitchwise direction and 85 computational points in the 
axial direction, with 49 computational points between leading and trailing 
edges. Figure 9 shows the computed constant velocity lines. The velocity is 
zero at the blade surface. The predicted Navier-Stokes flow field clearly 
shows the growth of the boundary layer over the blade surfaces. The growth on 
the rear suction side is particularly strong, but the adopted computational 
grid is not fine enough to allow the computation of recirculating flows. Fur- 
ther Navier-Stokes computations with higher grid refinement close to the blade 
rear suction side seem to be necessary before the final acceptance of the 
blade shape. 

It is to be noted that oscillations of very low amplitude continue inde- 
finitely, because of the wave reflections from the boundaries of the computa- 
tional domain, but the influence of these oscillations on the blade pressure 
distribution is negligible. 

Finally, the numerical method was applied to the computation of the rotor 
flow perturbed by the presence of downstream bodies. The computations were 
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performed by starting from Initial steady conditions and then superimposing to 
the steady value for the exit pressure a sinusoidal dependence of the form 

p 2 <t> = 98 000 [1 + 0.001 sin (0.001 t) ] N/m 2 
where t Is the time in sec. 

The successive Integration of the Euler equations yields (after a certain 
number of time steps) a periodic flow condition, where' the rotor pressure 
field follows an almost sinusoidal oscillation. Figure 10 shows an example of 
rotor pressure fluctuations on the blade suction and pressure side at the 
80-percent chord position. 

The computed pressure fluctuations are quite close to those expected, per- 
haps very close to those specified. at the outflow boundary. The method appears 
to be able to simulate sinusoidal 'pressure oscillations without introducing 
strong nonphysical dissipations. Obviously, a better test case is required to 
properly determine the prediction capability of the code. Such a comparison 
will be possible only when experimental data are available in enough 
detail. 


CONCLUSIONS 

This paper has presented a time-accurate numerical method that can prop- 
erly describe the steady and unsteady fan rotor flow fields in the 
blade-to-blade plane. The numerical integration algorithm reduces the computa- 
tional cost of an unsteady simulation involving high-frequency disturbances In 
both CPU time and memory requirements. During the computations performed on 
the NASA-Ames Cray Y-MP supercomputer, the advance of the flow equations in a 
computational grid made up of about 2000 grid points performing 10 000 time 
steps requires about 300 sec CPU time when solving the Nav1er-Stok.es equations 
and about 200 sec CPU time when solving the Euler equations. The memory 
required is less than 0.3 megawords. Thanks to the low computational cost and 
the good accuracy obtained during the simulations, the method appears to be a 
useful device to modify blade shape for better steady, aerodynamic perfor- 
mances. The code also Indicates some ways to Improve the unsteady aerodynamic 
performance . 
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APPENDIX - SYMBOLS 


A V artificial viscosity coefficient 

a speed of sound 

abs absolute 

C constant 

C e kinematic energy constant 

Cy eddy viscosity constant 

CFL Courant-Frledrlch-Levy stability criteria 

D wall term 

D diffusion vector 

dij Kronecker delta 

dr mesh size 

dt time step 

E wall term 

e kinetic energy dissipation rate 

F compressibility term 

f function 

f unknown vector 

G compressibility term 

h enthalpy 

K turbulence kinetic energy 

Le length scale 

M Mach number 

P turbulence kinetic energy production 

Pr Prandtl number 

p pressure 

q heat flux vector 
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Rt turbulence Reynolds number 
Sc Schmidt number 

T residual vector 

Tu turbulence level 

t time 

u velocity component 

V velocity 

(5 flow angle 

Y specific heat ratio 

ic heat transfer coefficient 

]i viscosity coefficient 

p density 

x viscous stress tensor 

Subscripts: 
t turbulent 

0 total 

1 inlet 

2 outlet 
Superscripts: 

m time level 

2 second order 

4 fourth order 
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FIGURE A . - STABLE TIME STEPS AT DIFFERENT COURANT-FR IEDR ICH-LEVY STABILITY CRITERION (CFL) NUMBERS. 
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FIGURE 7. - COMPARISON OF PRESCRIBED AND COMPUTED 
SURFACE MACH NUMBER DISTRIBUTIONS. 



FIGURE 8. - NAVIER-STOKES COMPUTATIONAL GRID. 
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